Continuous and discrete fields¶
Definition¶
In HySoP, a continuous field
(Field) is an abstract object
which represents the usual scalar, vector or tensor field, in
mathematical sense, i.e some function which associates a scalar, a
vector of any length tensor or a tensor of any order to each point of the space.
Such objects are used as input and/or output variables for continuous operators and must be defined with at least:
a
Domain: the physical domain of definition of the field,a name, compulsory (required for i/o).
a shape, optional, as a tuple of integers
Scalar fields are built by default. Use is_vector or nb_component optional arguments to build n-components fields (n=domain.dimension for vector fields).
>>> from hysop import Field, Box
>>> # build a domain
>>> box = Box(dim=3)
>>> # Scalar field
>>> scal = Field(domain=box, name='s1')
>>> # Vector field
>>> vec = Field(domain=box, name='v1', is_vector=True)
>>> # n-component field
>>> n = 4
>>> vec_n = Field(domain=box, name='vn', nb_components=n)
>>> # nb_component = domain.dimension is equivalent to is_vector=True
>>> to2 = Field(domain=box, name='v2', shape=(3,3))
Obviously one needs to associate one or more topologies (CartesianTopology, see MPI topologies and space discretisation) to such a field, to be able to apply operators or numerical methods. Remind that ‘topologies’ means either the definition of a global space resolution and of the way data are distributed through mpi processes
>>> from hysop.topology.cartesian_topology import CartesianTopology
>>> from hysop.tools.parameters import CartesianDiscretization
>>> d3d = CartesianDiscretization([33, ]*3, default_boundaries=True)
>>> topo = CartesianTopology(box, d3d)
>>> # build a discretisation of vec on d3d and distribute its data according to 'topo'
>>> vd = vec.discretize(topo)
Each time ‘discretize’ function is called, a new entry is added to the continuous field attribute named discreteFields, which maps a topology to a DiscreteField. Therefore, a discrete field will handle local arrays of values
of the field on the local mesh.
To summarize, a discrete field (DiscreteField) is the association of:
a topology (
CartesianTopology),a list of numpy arrays (one per component), values of the field on the local mesh.
For example
>>> d3d_with_ghosts = CartesianDiscretization([33, ]*3, ghosts=[2, ]*3, default_boundaries=True)
>>> topo_g = CartesianTopology(box, d3d_with_ghosts)
>>> # build a discretisation of vec on d3d and distribute its data according to 'topo'
>>> vd2 = vec.discretize(topo_g)
>>> for df in vd:
... print(df.resolution)
...
[33 33 33]
[33 33 33]
[33 33 33]
>>> for df in vd2:
... print(df.resolution)
...
[37 37 37]
[37 37 37]
[37 37 37]
will display the local resolution of each component of the discretisations ‘topo’ and ‘topo_g’ of vec. Their values will depend either on the global resolution, the number of mpi processes on which the program run and on the current process rank.
Notice that the exact behavior of
>>> vd = vec.discretize(topo)
is
check if a discretisation of vec on topo already exists. If so, returns this discretisation,
else discretizes vec on topo and returns the discrete field.
Remark:
Usually, in a ‘complete’ problem, the user does not have to bother with the topologies’ stuff and the fields discretisations. All this is hidden from user : one just has to define some operators and the fields they work with, how those fields are initialized (see below) , bricks of the global problem, associated with some required global resolutions and some solving methods. Then the setup process of the problem/operators will automatically performs the topologies creation and the discretisation of the fields. See operators/problem documentation for more details or the examples in HySoP/Examples.
How to initialize a field?¶
Fields values are computed with a user-defined function. To do so, continuous field has an ‘initialize’ method.
A user-defined function myfunction must be a function of time and space, filling n numpy arrays, n being the number of components of the field.
Therefore, to initialize vec, run
>>> def myfunction(data, coords, component):
... (x, y, z) = coords
... data[...] = x*y*z
...
>>> vd.initialize(myfunction)
One can add extra parameters to the function, for instance a time parameter
>>> def myfunction(data, coords, component, t=0.):
... (x, y, z) = coords
... data[...] = x*y*z*t
...
>>> vd.initialize(myfunction) # Use default value for parameter t
>>> vd.initialize(myfunction, t=0.1)
Remarks:
in 2d, z argument is removed, in 1d, only x remains (with t).
Res must be a list of numpy arrays of length equal to the number of component of the field to be initialized.
The number of returned values in case 1 must be equal to the number of component of the field to be initialized.
Case1 function takes scalar arguments as inputs and so will be vectorized during initialization (about vectorization in numpy, see http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.vectorize.html), while case2 function works on numpy arrays. Case 2 must be prefered when possible for performance reasons.
When vectorization is needed, vectorize_formula=True must be set during definition of the field.
Here is an example where a vector field ‘velo’ on a 2d domain is initialized with a vectorized function such that
and a vector field \(\omega\) on a 3d domain with a ‘numpy-ready’ function. Notice that this last case is what we call ‘Taylor-Green’ in many of our examples.
The corresponding implementation will be
>>> from numpy import cos, sin
>>> def func2d(data, coords, component, time=0.):
... (x, y) = coords
... if component == 0:
... data[...] = cos(x)
... if component == 1:
... data[...] = cos(y) + x*t
...
>>> # build a 2d domain
>>> box2d = Box(length=[1., 1.])
>>> # vector field
>>> v2 = Field(domain=box2d, name='v2', is_vector=True)
>>> d2d = CartesianDiscretization([33, ]*2, default_boundaries=True)
>>> topo2d = CartesianTopology(box2d, d2d)
>>> v2d = v2.discretize(topo2d)
>>> t=3.2
>>> v2d.initialize(func2d, time=t)
>>> def w_TG(data, coords, component, time):
... (x, y, z) = coords
... if component == 0:
... data[...] = - cos(x) * sin(y) * sin(z)
... if component == 1:
... data[...] = - sin(x) * cos(y) * sin(z)
... if component == 2:
... data[...] = 2. * sin(x) * sin(y) * cos(z)
...
>>> # build a 3d domain
>>> box3d = Box(length=[1., 1., 1.])
>>> # vector field
>>> v3 = Field(domain=box3d, name='v3', is_vector=True)
>>> d3d = CartesianDiscretization([33, ]*3, default_boundaries=True)
>>> topo3d = CartesianTopology(box3d, d3d)
>>> v3d = v3.discretize(topo3d)
>>> t=3.2
>>> v3d.initialize(w_TG, time=t)
Notice that the local mesh resolution or the data distribution are hidden from the user-defined function.
Other ways to initialize a field¶
Copy an existing field
>>> v = Field(domain=box2d, name='v', is_vector=True)
>>> vd = v.discretize(topo2d)
>>> vd.initialize(v2d.data) # initialise from v2d datas
Random initialization
>>> vd = vd.randomize()
Fields’ data access¶
Some examples on how to get access or set values of the fields
>>> d2dg = CartesianDiscretization([33, ]*2, ghosts=[2, ]*2, default_boundaries=True)
>>> v2d = v.discretize(CartesianTopology(box2d, d2dg))
>>> # Access to all discretization of a given continuous field
>>> for t,f in v[0].discrete_fields.items():
... print(type(t), f.data[0].shape)
...
<class 'hysop.topology.cartesian_topology.CartesianTopology'> (33, 33)
<class 'hysop.topology.cartesian_topology.CartesianTopology'> (37, 37)
>>> # Acces a given index of second component of a discrete field
>>> vd.data[1][2,3] = 3
>>> print(vd.data[1][2,3])
3.0
Writing/reading fields to/from file¶
Fields can be written to or read from hdf files thanks to numpy interface.
>>> import numpy as np
>>> file = "./dumped_example.npy"
>>> vd.data[0].dump(file)
>>> vd.data[1][...] = np.load(file)
>>> import os; os.remove(file) # Cleanning
h5py module can be used instead. For efficient IO, the usage of hysop’s IO operators is recommended (see I/O management: tools and conventions in HySoP).
Other useful functions¶
norm: \((\sum_{i,j,...} (|v_d(i,j,...)|^p) )^{1/p}\) (default is L2-norm)
Sum over all grid points of the values of component d of the discrete field. Computed with:
>>> vd.data[0][...] = 2.
>>> vd.data[1][...] = 1.
>>> print(vd.norm())
[2. 1.]
integrate: \(\int_\Omega v_d(i,j,...) d\boldsymbol{x}\) (default is L2-norm)
Integrate over domain, using sum over all grid points and scaled by elementary volume dx.
>>> print(vd.integrate())
[2. 1.]
